##############################################################################
##############################################################################
#
#				"Do Natural Disasters Help the Environment? How Voters 
#         Respond and What That Means"
#
#       Leonardo Baccini and Lucas Leemann
#
#       Contact: leemann@ipz.uzh.ch
#
#       Version: March 2020
#	 
##############################################################################
##############################################################################



##############################################################################
# Packages

library(foreign) # CRAN v0.8-72
library(maptools) # CRAN v0.9-9
library(RColorBrewer) # CRAN v1.1-2
library(Matching) # CRAN v4.9-6
library(ebal) # CRAN v0.1-6
library(xtable) # CRAN v1.8-4
library(survey) # CRAN v3.36
library(texreg) # CRAN v1.36.23
library(rgdal) # CRAN v1.4-8
library(dplyr) # [github::tidyverse/dplyr] v0.8.99.9002
library(ggplot2) # CRAN v3.3.0
library(raster) # CRAN v3.0-12
library(gpclib) # CRAN v1.5-6
library(rgeos) # CRAN v0.5-2


##############################################################################
##############################################################################
#
# Replication Code for Paper & Appendix
#
##############################################################################
##############################################################################


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Data manipulation and estimation of models for Table 1 in paper
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


data.st <- read.dta("Data/all data for CH.dta")
# take out the three environmental votes that are not related to climate change
data.st <- data.st[data.st$vote!=226 & data.st$vote!=227 & data.st$vote!=251,]
vdata <- data.st
vdata <- cbind(vdata$yes_vote_meanD,vdata$yes_p_per,vdata$treated, vdata$kuenstl,vdata$gras,vdata$gebuesch,vdata$baum,vdata$vlose,vdata$wasser,vdata$alti_mun_km,
               vdata$rainfall_surf,vdata$rainfall_1000,vdata$rainfall_1000_2,vdata$Punktflaeche, vdata$SPS_vshare, vdata$CVP_vshare, vdata$GPS_vshare, vdata$FDP_vshare, 
               vdata$SVP_vshare, vdata$vote, vdata$steep_percent, vdata$voteL1, vdata$voteL2, vdata$voteL3, vdata$voteL4, vdata$voteL7, vdata$voteL8, vdata$voteL9,
               vdata$voteL10,vdata$voteL11, vdata$damageLOW, vdata$damageMIDDLE, vdata$damageHIGH, vdata$elapsedTime, vdata$small_left_vshare, vdata$Gde_typen22, 
               vdata$Kantonsnummer, vdata$Area_50years_km2, vdata$gdenr, vdata$Tert_share, vdata$Tert_sekun_share, vdata$Steuerb_Einkommen, vdata$Steuerpflichtige)
               #vdata$turnouta)
#vdata <- na.omit(vdata)
vdata <- data.frame(vdata)
colnames(vdata) <- c("yes_vote_meanD","yes_p_per","treated", "kuenstl", "gras", "gebuesch", "baum", "vlose", "wasser", "alti_mun_km", "rainfall_surf", "rainfall_1000", 
                     "rainfall_1000_2", "Punktflaeche", "SPS_vshare", "CVP_vshare", "GPS_vshare", "FDP_vshare", "SVP_vshare", "vote", "steepness", "voteL1", "voteL2", 
                     "voteL3", "voteL4", "voteL7", "voteL8", "voteL9", "voteL10", "voteL11", "damage1", "damage2","damage3","elapsedTime", "small_left_vshare",
                     "Gde_typen22", "cantonnr", "Area_50years_km2", "gdenr", "Tertiar", "TertiarSekund", "income", "TaxPopulation")#, "Turnout")




# Change units on some of the geographical variables
vdata$kuenst.sh <- 100* vdata$kuenstl/vdata$Punktflaeche
vdata$wasser.sh <- 100* vdata$wasser/vdata$Punktflaeche
vdata$gras.sh <- 100* vdata$gras/vdata$Punktflaeche
vdata$gebuesch.sh <- 100* vdata$gebuesch/vdata$Punktflaeche
vdata$baum.sh <- 100* vdata$baum/vdata$Punktflaeche
vdata$vlose.sh <- 100* vdata$vlose/vdata$Punktflaeche
vdata$alti_mun_m <- 1000* vdata$alti_mun_km
vdata$alti_mun_m2 <- vdata$alti_mun_m^2
vdata$rainfall_surf_1000 <- 1000*vdata$rainfall_surf
vdata$steep_per <- 100 * vdata$steepness


mod1 <- lm(yes_p_per ~  treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare + SVP_vshare +factor(vote) + factor(gdenr),
           data=vdata)
mod2 <- lm(yes_p_per ~ treated +  rainfall_1000  + vlose.sh + wasser.sh + gras.sh + 
             kuenst.sh  +factor(vote) + factor(gdenr), data=vdata)
mod3 <- lm(yes_p_per ~ treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare+ SVP_vshare+  rainfall_1000  +
             vlose.sh + wasser.sh + gras.sh + kuenst.sh  + factor(vote)+ factor(gdenr), data=vdata)

screenreg(list(mod1,mod2,mod3), stars=c(0.01,0.05,0.1), omit="factor", reorder.coef=c(2:7,8:12,1), 
          custom.coef.names = c("Intercept", "Flooded","Green Party %","Social Democrats %" , "Christian Democrats %", 
                                "Liberal Democrats %", "Swiss People's Party %","Rainfall", "No vegetation",
                                "Share of Water", "Share of Gras", "Artificial"), 
          digits=2) 

sink("Replication of Tables A.txt")
print("######")
print("######")
print("Table 1")
print("######")
print("######")
screenreg(list(mod1,mod2,mod3), stars=c(0.01,0.05,0.1), omit="factor", reorder.coef=c(2:7,8:12,1), 
          custom.coef.names = c("Intercept", "Flooded","Green Party %","Social Democrats %" , "Christian Democrats %", 
                                "Liberal Democrats %", "Swiss People's Party %","Rainfall", "No vegetation",
                                "Share of Water", "Share of Gras", "Artificial"), 
          digits=2) 

sink()






# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Run models for Table 2 in paper
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

vdata1 <- na.omit(vdata)
treatment   <- vdata1$treated 
X           <- cbind(vdata1$vlose.sh, vdata1$wasser.sh, vdata1$gras.sh, vdata1$kuenst.sh, vdata1$alti_mun_m, vdata1$rainfall_surf_1000, vdata1$steepness, 
                     vdata1$SPS_vshare ,vdata1$CVP_vshare,vdata1$GPS_vshare,vdata1$FDP_vshare,vdata1$SVP_vshare,vdata1$voteL1, vdata1$voteL2, vdata1$voteL3, 
                     vdata1$voteL4, vdata1$voteL8, vdata1$voteL9, vdata1$voteL10, vdata1$Area_50years_km2)
colnames(X) <- c("no vegetation" ,"wasser", "gras","kuenstl",  "altitude", "rainfall", "steepness", "SPS", "CVP", "GPS", "FDP", "SVP", "vote1", "vote2", 
                 "vote3", "vote4", "vote8", "vote9", "vote10", "Area_50years_km2")

# Entropy Balancing
eb.out <- ebalance(Treatment=treatment, X=X)
weight.L <- rep(1,dim(X)[1])
weight.L[treatment==0] <- eb.out$w

Yl  <- vdata1$yes_p_per
mb.ent.before <- MatchBalance(treatment~ vdata1$Area_50years_km2+vdata1$vlose.sh+ vdata1$wasser.sh+ vdata1$gras.sh  + vdata1$kuenst.sh + vdata1$alti_mun_m + 
                                vdata1$rainfall_surf_1000 + vdata1$steepness + vdata1$SPS_vshare + vdata1$CVP_vshare + vdata1$GPS_vshare + vdata1$FDP_vshare + 
                                vdata1$SVP_vshare, weights=rep(1,length(weight.L)), ks = FALSE)
mb.ent.after <- MatchBalance(treatment~ vdata1$Area_50years_km2+vdata1$vlose.sh+ vdata1$wasser.sh+ vdata1$gras.sh  + vdata1$kuenst.sh  + vdata1$alti_mun_m + vdata1$rainfall_surf_1000 + 
                               vdata1$steepness + vdata1$SPS_vshare + vdata1$CVP_vshare + vdata1$GPS_vshare + vdata1$FDP_vshare + vdata1$SVP_vshare, 
                             weights=weight.L, ks = FALSE)
bloc.pre <- baltest.collect(matchbal.out=mb.ent.before, var.names=c( "Flooding Risk","Surface: % No vegetation", "Surface: % Water","Surface: % Gras","Surface: % Artificial","Altitude (in m)", "Rainfall (per sqkm)", "Steepness in %", "Social Democrats (%)", "Christian Democrats (%)", "Greens (%)", "Liberals (%)", "Swiss People's Party (%)"), after=FALSE)
bloc.post <- baltest.collect(matchbal.out=mb.ent.after, var.names=c("Flooding Risk","Surface: % No vegetation", "Surface: % Water", "Surface: % Gras","Surface: % Artificial","Altitude (in m)", "Rainfall (per sqkm)", "Steepness in %", "Social Democrats (%)", "Christian Democrats (%)", "Greens (%)", "Liberals (%)", "Swiss People's Party (%)"), after=FALSE)
rowN <- rownames(bloc.post)
blocL <- matrix(NA,length(rowN),6)
rownames(blocL) <- rowN
blocL[,c(1:3)] <- round(bloc.pre[,c(1,2,6)],2)
blocL[,c(4:5)] <- round(bloc.post[,c(1,2)],2)


bloc.post
# get treatment effect estimate
dat <- data.frame(Yl,vdata1$cantonnr, weight.L,treatment,vdata1$SPS_vshare,vdata1$GPS_vshare,vdata1$CVP_vshare,vdata1$FDP_vshare,vdata1$SVP_vshare,vdata1$elapsedTime,vdata1$gdenr, vdata1$Tertiar, vdata1$TertiarSekund, vdata1$income)
des1 <- svydesign(id=~1,weights=~weight.L, data= dat)
mod1 <- svyglm(Yl ~ treatment , design = des1)

summary(mod1)
blocL
for (i in 1:dim(blocL)[1]){
  # bloc.post[,6] holds the t-values and they are all essentially 1, i.e. showing no difference
  if (round(bloc.post[i,6],2)==1.00) blocL[i,6] <- "\\checkmark"
}

sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("### Table 2 (ATT)")
print("######")
print("######")
summary(mod1)
blocL
sink()




# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Run models for Table 3 in paper
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

vdata1$left <- vdata1$SPS_vshare + vdata1$GPS_vshare + vdata1$small_left_vshare 
vdata1$leftINTtreat <- vdata1$left * vdata1$treated
vdata1$elapsedTime[treatment==0] <- 0
vdata1$timeINTtreat <- vdata1$elapsedTime * vdata1$treated
dat <- data.frame(Yl,vdata1$cantonnr, weight.L,treatment,vdata1$leftINTtreat,vdata1$left,
                  vdata1$SPS_vshare,vdata1$GPS_vshare,vdata1$CVP_vshare,vdata1$FDP_vshare,
                  vdata1$SVP_vshare,vdata1$elapsedTime,vdata1$timeINTtreat, vdata1$gdenr, 
                  vdata1$Tertiar, vdata1$TertiarSekund, vdata1$income)
des1 <- svydesign(id=~1,weights=~weight.L, data= dat)
mod1 <- svyglm(Yl ~ treatment , design = des1)
mod2 <- svyglm(Yl ~ treatment + vdata1.Tertiar, design = des1)
mod3 <- svyglm(Yl ~ treatment + vdata1.Tertiar+ vdata1.Tertiar:treatment, design = des1)


screenreg(list(mod1, mod2,mod3),  stars=c(0.01,0.05,0.1))
sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table 3 (Heterogeneity in Space)")
print("######")
print("######")
screenreg(list(mod1, mod2,mod3),  stars=c(0.01,0.05,0.1))
sink()


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Run models for Table 4 in paper
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

mod4 <- svyglm(Yl ~ treatment , design = des1)
mod5 <- svyglm(Yl ~ treatment +  vdata1.elapsedTime + vdata1.timeINTtreat, design = des1)
screenreg(list(mod4,mod5),  stars=c(0.01,0.05,0.1), custom.coef.names = c("Constant","Treatment","Time betw. Flood and Vote"), custom.model.names = c("Model V","Model VIII"))

sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table 4 (Heterogeneity in Time)")
print("######")
print("######")
screenreg(list(mod4,mod5),  stars=c(0.01,0.05,0.1), custom.coef.names = c("Constant","Treatment","Time betw. Flood and Vote"), custom.model.names = c("Model V","Model VIII"))
sink()



# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Figure 2
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Generating Quantities, then plotting them 
# Education
vL <- vcov(mod3)
bL <- coef(mod3)
BETA3 <- mvrnorm(n = 5000,mu=bL, Sigma = vL)
head(BETA3)
educ <- seq(1,100,by = 1)/100
Yhat <- BETA3[,c(2,4)] %*% t(cbind(1,educ))

Y.3 <- matrix(NA,100,3)
for (i in 1:100){
  Y.3[i,] <- quantile(Yhat[,i], c(0.025, 0.5, 0.975))
}
# time
vL <- vcov(mod5)
bL <- coef(mod5)
BETA7 <- mvrnorm(n = 5000,mu=bL, Sigma = vL)
head(BETA7)
timeL <- seq(1,120,by = 1)/12
Yhatt <- BETA7[,c(2,3)] %*% t(cbind(1,timeL))

Y.7 <- matrix(NA,120,3)
for (i in 1:120){
  Y.7[i,] <- quantile(Yhatt[,i], c(0.025, 0.5, 0.975))
}
# Plot
pdf("Marginaleffects.pdf",width=10,height=5)
par(mfrow=c(1,2))
plot(educ,colMeans(Yhat), col="blue", type = "l",
     lwd=3, bty="n", ylab="Estimated Effect of Exposure",
     xlab="Share of People with Tertiary Education")
points(educ,Y.3[,1], type="l", lty=3,
       col="blue", lwd=2) 
points(educ,Y.3[,3], type="l", lty=3,
       col="blue", lwd=2) 

for (i in 1:5000){
  points(educ,Yhat[i,], type="l",
         col=rgb(0,0,255,2,maxColorValue = 255), lwd=2) 
}
# 2nd plot
plot(timeL,colMeans(Yhatt), col="blue", type = "l",
     lwd=3, bty="n", ylab="Estimated Effect of Exposure",
     xlab="Months since Exposure")
points(timeL,Y.7[,1], type="l", lty=3,
       col="blue", lwd=2) 
points(timeL,Y.7[,3], type="l", lty=3,
       col="blue", lwd=2) 

for (i in 1:5000){
  points(timeL,Yhatt[i,], type="l",
         col=rgb(0,0,255,2,maxColorValue = 255), lwd=2) 
}
dev.off()


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Significance test in FN 12
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


set.seed(123)
# significance check:
vL <- vcov(mod5)
bL <- coef(mod5)
BETA <- mvrnorm(n = 10000,mu=bL, Sigma = vL)
# After 10 months no significant difference
x0 <- c(1,0,0)
x1 <- c(1,1,10)
p0 <- BETA %*% x0
p1 <- BETA %*% x1
diff01 <- p0 - p1

sort(diff01)[c(500, 9501)] # 90
sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Footnote 12, 90% CI")
print("######")
print("######")
sort(diff01)[c(500, 9501)] # 90
sink()


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Table A2
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table A1")
print("######")
print("######")
print("See above where Table 1 is replicated - they are identical.")
sink()


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Table A2
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

data1 <- read.csv("Data/pnl_leo_w4.csv")
data1$cc.worry <- rep(NA,dim(data1)[1])
data1$Q16n <- factor(data1$Q16, levels=c("Gar nicht besorgt",
                                         "Eher nicht besorgt",
                                         "Eher besorgt",
                                         "Sehr besorgt",
                                         "Weiss nicht / keine Antwort"))

data1$Q2n <- factor(data1$Q2, levels=c("Gar nicht besorgt",
                                       "Eher nicht besorgt",
                                       "Eher besorgt",
                                       "Sehr besorgt",
                                       "Weiss nicht / keine Antwort"))
data1$cc.worry[is.na(data1$Q2n)] <- data1$Q16n[is.na(data1$Q2n)]
data1$cc.worry[is.na(data1$Q16n)] <- data1$Q2n[is.na(data1$Q16n)]

pdf("Surveyresponse.pdf",width=10,height=15)
par(mar=c(3,15,1,1))
barplot(table(data1$Q4),horiz = TRUE, las=2,
        names.arg = c("Deregulation of Construction Rules",
                      "Never wondered about this",
                      "No clear `reason`",
                      "Ineffective Bureaucracy",
                      "Climate Change",
                      "Don't know"),
        col="dodgerblue3",
        border=NA)
dev.off()


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Table A2
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

data.st <- read.dta("Data/all data for CH.dta")
# take out the three environmental votes that are not related to climate change
data.st <- data.st[data.st$vote!=226 & data.st$vote!=227 & data.st$vote!=251,]
vdata <- data.st
vdata <- cbind(vdata$yes_vote_meanD,vdata$yes_p_per,vdata$treated, vdata$kuenstl,vdata$gras,vdata$gebuesch,vdata$baum,vdata$vlose,vdata$wasser,vdata$alti_mun_km,
               vdata$rainfall_surf,vdata$rainfall_1000,vdata$rainfall_1000_2,vdata$Punktflaeche, vdata$SPS_vshare, vdata$CVP_vshare, vdata$GPS_vshare, vdata$FDP_vshare, 
               vdata$SVP_vshare, vdata$vote, vdata$steep_percent, vdata$voteL1, vdata$voteL2, vdata$voteL3, vdata$voteL4, vdata$voteL7, vdata$voteL8, vdata$voteL9,
               vdata$voteL10,vdata$voteL11, vdata$damageLOW, vdata$damageMIDDLE, vdata$damageHIGH, vdata$elapsedTime, vdata$small_left_vshare, vdata$Gde_typen22, 
               vdata$Kantonsnummer, vdata$Area_50years_km2, vdata$gdenr, vdata$Tert_share, vdata$Tert_sekun_share, vdata$Steuerb_Einkommen, vdata$Steuerpflichtige,vdata$turnouta)
#vdata <- na.omit(vdata)
vdata <- data.frame(vdata)
colnames(vdata) <- c("yes_vote_meanD","yes_p_per","treated", "kuenstl", "gras", "gebuesch", "baum", "vlose", "wasser", "alti_mun_km", "rainfall_surf", "rainfall_1000", 
                     "rainfall_1000_2", "Punktflaeche", "SPS_vshare", "CVP_vshare", "GPS_vshare", "FDP_vshare", "SVP_vshare", "vote", "steepness", "voteL1", "voteL2", 
                     "voteL3", "voteL4", "voteL7", "voteL8", "voteL9", "voteL10", "voteL11", "damage1", "damage2","damage3","elapsedTime", "small_left_vshare",
                     "Gde_typen22", "cantonnr", "Area_50years_km2", "gdenr", "Tertiar", "TertiarSekund", "income", "TaxPopulation", "Turnout")




# Change units on some of the geographical variables
vdata$kuenst.sh <- 100* vdata$kuenstl/vdata$Punktflaeche
vdata$wasser.sh <- 100* vdata$wasser/vdata$Punktflaeche
vdata$gras.sh <- 100* vdata$gras/vdata$Punktflaeche
vdata$gebuesch.sh <- 100* vdata$gebuesch/vdata$Punktflaeche
vdata$baum.sh <- 100* vdata$baum/vdata$Punktflaeche
vdata$vlose.sh <- 100* vdata$vlose/vdata$Punktflaeche
vdata$alti_mun_m <- 1000* vdata$alti_mun_km
vdata$alti_mun_m2 <- vdata$alti_mun_m^2
vdata$rainfall_surf_1000 <- 1000*vdata$rainfall_surf
vdata$steep_per <- 100 * vdata$steepness


mod1 <- lm(Turnout ~  treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare + SVP_vshare +factor(vote) + factor(gdenr),
           data=vdata)
mod2 <- lm(Turnout ~ treated +  rainfall_1000 + vlose.sh + wasser.sh + gras.sh + 
             kuenst.sh  +factor(vote) + factor(gdenr), data=vdata)
mod3 <- lm(Turnout ~ treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare+ SVP_vshare+  rainfall_1000 + 
             vlose.sh + wasser.sh + gras.sh + kuenst.sh  + factor(vote) + factor(gdenr), data=vdata)

screenreg(list(mod1,mod2,mod3), stars=c(0.01,0.05,0.1), omit="factor", reorder.coef=c(2:7,8:12,1),#c(2:10,16,11:15,1), 
       custom.coef.names = c("Intercept", "Flooded","Green Party %","Social Democrats %" , "Christian Democrats %", 
                             "Liberal Democrats %", "Swiss People's Party %","Rainfall", "No vegetation","Share of Water", 
                             "Share of Gras", "Artificial"), 
       digits=2) 

sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table A2")
print("######")
print("######")
screenreg(list(mod1,mod2,mod3), stars=c(0.01,0.05,0.1), omit="factor", reorder.coef=c(2:7,8:12,1),#c(2:10,16,11:15,1), 
          custom.coef.names = c("Intercept", "Flooded","Green Party %","Social Democrats %" , "Christian Democrats %", 
                                "Liberal Democrats %", "Swiss People's Party %","Rainfall", "No vegetation","Share of Water", 
                                "Share of Gras", "Artificial"), 
          digits=2) 
sink()




# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Table A3
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

data.st <- read.dta("Data/all data for CH ROBUSTNESS FORWARD1.dta")
data.st <- data.st[data.st$vote!=226 & data.st$vote!=227 & data.st$vote!=251,]
vdata <- data.st
vdata <- cbind(vdata$yes_vote_meanD,vdata$yes_p_per,vdata$treated,vdata$treated_after, vdata$kuenstl,vdata$gras,vdata$gebuesch,vdata$baum,vdata$vlose,vdata$wasser,vdata$alti_mun_km,
               vdata$rainfall_surf,vdata$rainfall_1000,vdata$rainfall_1000_2,vdata$Punktflaeche, vdata$SPS_vshare, vdata$CVP_vshare, vdata$GPS_vshare, vdata$FDP_vshare, 
               vdata$SVP_vshare, vdata$vote, vdata$steep_percent, vdata$voteL1, vdata$voteL2, vdata$voteL3, vdata$voteL4, vdata$voteL7, vdata$voteL8, vdata$voteL9,
               vdata$voteL10,vdata$voteL11, vdata$damageLOW, vdata$damageMIDDLE, vdata$damageHIGH, vdata$elapsedTime, vdata$small_left_vshare, vdata$Gde_typen22, 
               vdata$Kantonsnummer, vdata$Area_50years_km2, vdata$gdenr, vdata$Tert_share, vdata$Tert_sekun_share, vdata$Steuerb_Einkommen, vdata$Steuerpflichtige)
vdata <- data.frame(vdata)
colnames(vdata) <- c("yes_vote_meanD","yes_p_per","treated", "treated_after", "kuenstl", "gras", "gebuesch", "baum", "vlose", "wasser", "alti_mun_km", "rainfall_surf", "rainfall_1000", 
                     "rainfall_1000_2", "Punktflaeche", "SPS_vshare", "CVP_vshare", "GPS_vshare", "FDP_vshare", "SVP_vshare", "vote", "steepness", "voteL1", "voteL2", 
                     "voteL3", "voteL4", "voteL7", "voteL8", "voteL9", "voteL10", "voteL11", "damage1", "damage2","damage3","elapsedTime", "small_left_vshare",
                     "Gde_typen22", "cantonnr", "Area_50years_km2", "gdenr", "Tertiar", "TertiarSekund", "income", "TaxPopulation")


vdata$kuenst.sh <- 100* vdata$kuenstl/vdata$Punktflaeche
vdata$wasser.sh <- 100* vdata$wasser/vdata$Punktflaeche
vdata$gras.sh <- 100* vdata$gras/vdata$Punktflaeche
vdata$gebuesch.sh <- 100* vdata$gebuesch/vdata$Punktflaeche
vdata$baum.sh <- 100* vdata$baum/vdata$Punktflaeche
vdata$vlose.sh <- 100* vdata$vlose/vdata$Punktflaeche
vdata$alti_mun_m <- 1000* vdata$alti_mun_km
vdata$alti_mun_m2 <- vdata$alti_mun_m^2
vdata$rainfall_surf_1000 <- 1000*vdata$rainfall_surf
vdata$steep_per <- 100 * vdata$steepness



mod01 <- lm(yes_p_per ~ treated_after + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare + SVP_vshare +factor(vote) + factor(gdenr),
            data=vdata)
mod1 <- lm(yes_p_per ~ treated + treated_after + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare + SVP_vshare +factor(vote) + factor(gdenr),
           data=vdata)
mod02 <- lm(yes_p_per ~  treated_after +  rainfall_1000 +  vlose.sh + wasser.sh + gras.sh + 
              kuenst.sh  +factor(vote) + factor(gdenr), data=vdata)

mod2 <- lm(yes_p_per ~ treated + treated_after +  rainfall_1000 + vlose.sh + wasser.sh + gras.sh + 
             kuenst.sh  +factor(vote)+ factor(gdenr), data=vdata)
mod03 <- lm(yes_p_per ~  treated_after + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare+ SVP_vshare+  rainfall_1000 + 
              vlose.sh + wasser.sh + gras.sh + kuenst.sh  + factor(vote) + factor(gdenr), data=vdata)
mod3 <- lm(yes_p_per ~ treated + treated_after + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare+ SVP_vshare+  rainfall_1000 +
             vlose.sh + wasser.sh + gras.sh + kuenst.sh  + factor(vote)+ factor(gdenr), data=vdata)


screenreg(list(mod01,mod1,mod02,mod2,mod03,mod3), omit="factor|(Intercept)",  reorder.coef=c(7,1,2:6,8,9:12), 
       custom.coef.names = c("Flooded After Vote","Green Party %","Social Democrats %" , "Christian Democrats %", 
                             "Liberal Democrats %", "Swiss People's Party %", "Flooded","Rainfall", "No vegetation",
                             "Share of Water", "Share of Gras", "Artificial"), digits=2,stars=c(0.01,0.05,0.1)) 

sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table A3")
print("######")
print("######")
screenreg(list(mod01,mod1,mod02,mod2,mod03,mod3), omit="factor|(Intercept)",  reorder.coef=c(7,1,2:6,8,9:12), 
          custom.coef.names = c("Flooded After Vote","Green Party %","Social Democrats %" , "Christian Democrats %", 
                                "Liberal Democrats %", "Swiss People's Party %", "Flooded","Rainfall", "No vegetation",
                                "Share of Water", "Share of Gras", "Artificial"), digits=2,stars=c(0.01,0.05,0.1)) 
sink()


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Table A4
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

vdata1 <- na.omit(vdata)
treatment   <- vdata1$treated_after # vlose.sh + wasser.sh + gras.sh + kuenst.sh
X           <- cbind(vdata1$vlose.sh, vdata1$wasser.sh, vdata1$gras.sh, vdata1$kuenst.sh, vdata1$alti_mun_m, vdata1$rainfall_surf_1000, vdata1$steepness, 
                     vdata1$SPS_vshare ,vdata1$CVP_vshare,vdata1$GPS_vshare,vdata1$FDP_vshare,vdata1$SVP_vshare,vdata1$voteL1, vdata1$voteL2, vdata1$voteL3, 
                     vdata1$voteL4, vdata1$voteL8, vdata1$voteL9, vdata1$voteL10, vdata1$Area_50years_km2)
colnames(X) <- c("no vegetation" ,"wasser", "gras","kuenstl",  "altitude", "rainfall", "steepness", "SPS", "CVP", "GPS", "FDP", "SVP", "vote1", "vote2", 
                 "vote3", "vote4", "vote8", "vote9", "vote10", "Area_50years_km2")

eb.out <- ebalance(Treatment=treatment, X=X)
weight.L <- rep(1,dim(X)[1])
weight.L[treatment==0] <- eb.out$w
Yl  <- vdata1$yes_p_per
mb.ent.before <- MatchBalance(treatment~ vdata1$Area_50years_km2+vdata1$vlose.sh+ vdata1$wasser.sh+ vdata1$gras.sh  + vdata1$kuenst.sh + vdata1$alti_mun_m + 
                                vdata1$rainfall_surf_1000 + vdata1$steepness + vdata1$SPS_vshare + vdata1$CVP_vshare + vdata1$GPS_vshare + vdata1$FDP_vshare + 
                                vdata1$SVP_vshare, weights=rep(1,length(weight.L)), ks = FALSE)
mb.ent.after <- MatchBalance(treatment~ vdata1$Area_50years_km2+vdata1$vlose.sh+ vdata1$wasser.sh+ vdata1$gras.sh  + vdata1$kuenst.sh  + vdata1$alti_mun_m + vdata1$rainfall_surf_1000 + 
                               vdata1$steepness + vdata1$SPS_vshare + vdata1$CVP_vshare + vdata1$GPS_vshare + vdata1$FDP_vshare + vdata1$SVP_vshare, 
                             weights=weight.L, ks = FALSE)
bloc.pre <- baltest.collect(matchbal.out=mb.ent.before, 
                            var.names=c( "Flooding Risk","Surface: % No vegetation", "Surface: % Water","Surface: % Gras",
                                         "Surface: % Artificial","Altitude (in m)", "Rainfall (per sqkm)", "Steepness in %", 
                                         "Social Democrats (%)", "Christian Democrats (%)", "Greens (%)", "Liberals (%)", 
                                         "Swiss People's Party (%)"), 
                            after=FALSE)
bloc.post <- baltest.collect(matchbal.out=mb.ent.after, var.names=c("Flooding Risk","Surface: % No vegetation", 
                                                                    "Surface: % Water", "Surface: % Gras",
                                                                    "Surface: % Artificial","Altitude (in m)", 
                                                                    "Rainfall (per sqkm)", "Steepness in %", "Social Democrats (%)", 
                                                                    "Christian Democrats (%)", "Greens (%)", "Liberals (%)", 
                                                                    "Swiss People's Party (%)"), 
                             after=FALSE)
rowN <- rownames(bloc.post)
blocL <- matrix(NA,length(rowN),6)
rownames(blocL) <- rowN
blocL[,c(1:3)] <- round(bloc.pre[,c(1,2,6)],2)
blocL[,c(4:5)] <- round(bloc.post[,c(1,2)],2)

dat <- data.frame(Yl,vdata1$cantonnr, weight.L,treatment,vdata1$SPS_vshare,
                  vdata1$GPS_vshare,vdata1$CVP_vshare,vdata1$FDP_vshare,
                  vdata1$SVP_vshare,vdata1$elapsedTime,vdata1$gdenr, vdata1$Tertiar, 
                  vdata1$TertiarSekund, vdata1$income)
des1 <- svydesign(id=~1,weights=~weight.L, data= dat)
mod1 <- svyglm(Yl ~ treatment , design = des1)



#xtable(blocL[,-6])
#summary(mod1)

for (i in 1:dim(blocL)[1]){
  if (round(bloc.post[i,6],2)==1.00) blocL[i,6] <- "\\checkmark"
}

sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table A4")
print("######")
print("######")
summary(mod1)
blocL
sink()


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Table A5
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

data.st <- read.dta("Data/all data for CH.dta")
data.st <- data.st[data.st$vote==226 | data.st$vote==227 | data.st$vote==251,]
vdata <- data.st
vdata <- cbind(vdata$yes_vote_meanD,vdata$yes_p_per,vdata$treated, vdata$kuenstl,vdata$gras,vdata$gebuesch,
               vdata$baum,vdata$vlose,vdata$wasser,vdata$alti_mun_km,vdata$rainfall_surf,vdata$rainfall_1000,
               vdata$rainfall_1000_2,vdata$Punktflaeche, vdata$SPS_vshare, vdata$CVP_vshare, vdata$GPS_vshare, 
               vdata$FDP_vshare, vdata$SVP_vshare, vdata$vote, vdata$steep_percent, vdata$voteL1, vdata$voteL2, 
               vdata$voteL3, vdata$voteL4,vdata$voteL5,vdata$voteL6, vdata$voteL7, vdata$voteL8, vdata$voteL9,
               vdata$voteL10,vdata$voteL11, vdata$damageLOW, vdata$damageMIDDLE, vdata$damageHIGH, vdata$elapsedTime, 
               vdata$small_left_vshare, vdata$Gde_typen22, vdata$Kantonsnummer, vdata$Area_50years_km2, vdata$gdenr, 
               vdata$Tert_sekun_share, vdata$Steuerb_Einkommen)
vdata <- data.frame(vdata)
colnames(vdata) <- c("yes_vote_meanD","yes_p_per","treated", "kuenstl", "gras", "gebuesch", "baum", "vlose", "wasser", 
                     "alti_mun_km", "rainfall_surf", "rainfall_1000", "rainfall_1000_2", "Punktflaeche", "SPS_vshare", 
                     "CVP_vshare", "GPS_vshare", "FDP_vshare", "SVP_vshare", "vote", "steepness", "voteL1", "voteL2", 
                     "voteL3", "voteL4",  "voteL5",  "voteL6", "voteL7", "voteL8", "voteL9", "voteL10", "voteL11", 
                     "damage1", "damage2","damage3","elapsedTime", "small_left_vshare", "Gde_typen22", "cantonnr", 
                     "Area_50years_km2","gdenr","Tertiar", "income")

vdata$kuenst.sh <- 100* vdata$kuenstl/vdata$Punktflaeche
vdata$wasser.sh <- 100* vdata$wasser/vdata$Punktflaeche
vdata$gras.sh <- 100* vdata$gras/vdata$Punktflaeche
vdata$gebuesch.sh <- 100* vdata$gebuesch/vdata$Punktflaeche
vdata$baum.sh <- 100* vdata$baum/vdata$Punktflaeche
vdata$vlose.sh <- 100* vdata$vlose/vdata$Punktflaeche
vdata$alti_mun_m <- 1000* vdata$alti_mun_km
vdata$alti_mun_m2 <- vdata$alti_mun_m^2
vdata$rainfall_surf_1000 <- 1000*vdata$rainfall_surf
vdata$steep_per <- 100 * vdata$steepness



mod1 <- lm(yes_p_per ~  treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare + SVP_vshare +factor(vote) + factor(gdenr),
           data=vdata)
mod2 <- lm(yes_p_per ~ treated +  rainfall_1000  + vlose.sh + wasser.sh + gras.sh + 
             kuenst.sh  +factor(vote) + factor(gdenr), data=vdata)
mod3 <- lm(yes_p_per ~ treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare+ SVP_vshare+  rainfall_1000 + 
             vlose.sh + wasser.sh + gras.sh + kuenst.sh  + factor(vote) + factor(gdenr), data=vdata)

screenreg(list(mod1,mod2,mod3), stars=c(0.01,0.05,0.1), omit="factor|(Intercept)", 
       groups=list("\\sc{Exposure}"=1,"\\sc{Vote Shares}"=2:6, "\\sc{Precipitation}"=7,  "\\sc{Surface}"=8:11), 
       reorder.coef=c(1:6,7:11), custom.coef.names = c("Flooded","Green Party %","Social Democrats %" , "Christian Democrats %", "Liberal Democrats %", 
                                                         "Swiss People's Party %","Rainfall", "No vegetation","Share of Water", 
                                                         "Share of Gras", "Artificial"), digits=2) 

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Table A6
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table A5")
print("######")
print("######")
screenreg(list(mod1,mod2,mod3), stars=c(0.01,0.05,0.1), omit="factor|(Intercept)", 
          groups=list("\\sc{Exposure}"=1,"\\sc{Vote Shares}"=2:6, "\\sc{Precipitation}"=7,  "\\sc{Surface}"=8:11), 
          reorder.coef=c(1:6,7:11), custom.coef.names = c("Flooded","Green Party %","Social Democrats %" , "Christian Democrats %", "Liberal Democrats %", 
                                                          "Swiss People's Party %","Rainfall", "No vegetation","Share of Water", 
                                                          "Share of Gras", "Artificial"), digits=2) 
sink()

vdata1 <- na.omit(vdata)
treatment   <- vdata1$treated
X           <- cbind(vdata1$vlose.sh, vdata1$wasser.sh, vdata1$gras.sh, vdata1$kuenst.sh, vdata1$alti_mun_m, 
                     vdata1$rainfall_surf_1000, vdata1$steepness, vdata1$SPS_vshare ,vdata1$CVP_vshare,
                     vdata1$GPS_vshare,vdata1$FDP_vshare,vdata1$SVP_vshare, vdata1$Area_50years_km2)
colnames(X) <- c("no vegetation" ,"wasser", "gras","kuenstl",  "altitude", "rainfall", "steepness", "SPS", 
                 "CVP", "GPS", "FDP", "SVP",  "Area_50years_km2")

eb.out <- ebalance(Treatment=treatment, X=X)
weight.L <- rep(1,dim(X)[1])
weight.L[treatment==0] <- eb.out$w
Yl  <- vdata1$yes_p_per
mb.ent.before <- MatchBalance(treatment~ vdata1$Area_50years_km2+vdata1$vlose.sh+ vdata1$wasser.sh+ vdata1$gras.sh  + 
                                vdata1$kuenst.sh + vdata1$alti_mun_m + vdata1$rainfall_surf_1000 + vdata1$steepness + 
                                vdata1$SPS_vshare + vdata1$CVP_vshare + vdata1$GPS_vshare + vdata1$FDP_vshare + 
                                vdata1$SVP_vshare, weights=rep(1,length(weight.L)), 
                              ks = FALSE)
mb.ent.after <- MatchBalance(treatment~ vdata1$Area_50years_km2+vdata1$vlose.sh+ vdata1$wasser.sh+ vdata1$gras.sh  + 
                               vdata1$kuenst.sh  + vdata1$alti_mun_m + vdata1$rainfall_surf_1000 + vdata1$steepness + 
                               vdata1$SPS_vshare + vdata1$CVP_vshare + vdata1$GPS_vshare + vdata1$FDP_vshare + vdata1$SVP_vshare, 
                             weights=weight.L, ks = FALSE)
bloc.pre <- baltest.collect(matchbal.out=mb.ent.before, var.names=c( "Flooding Risk","Surface: % No vegetation", 
                                                                     "Surface: % Water","Surface: % Gras","Surface: % Artificial",
                                                                     "Altitude (in m)", "Rainfall (per sqkm)", "Steepness in %", 
                                                                     "Social Democrats (%)", "Christian Democrats (%)", "Greens (%)", 
                                                                     "Liberals (%)", "Swiss People's Party (%)"), 
                            after=FALSE)
bloc.post <- baltest.collect(matchbal.out=mb.ent.after, var.names=c("Flooding Risk","Surface: % No vegetation", "Surface: % Water", 
                                                                    "Surface: % Gras","Surface: % Artificial","Altitude (in m)", 
                                                                    "Rainfall (per sqkm)", "Steepness in %", "Social Democrats (%)", 
                                                                    "Christian Democrats (%)", "Greens (%)", "Liberals (%)", 
                                                                    "Swiss People's Party (%)"), 
                             after=FALSE)
rowN <- rownames(bloc.post)
blocL <- matrix(NA,length(rowN),6)
rownames(blocL) <- rowN
blocL[,c(1:3)] <- round(bloc.pre[,c(1,2,6)],2)
blocL[,c(4:5)] <- round(bloc.post[,c(1,2)],2)
dat <- data.frame(Yl,vdata1$cantonnr, weight.L,treatment)
des1 <- svydesign(id=~1,weights=~weight.L, data= dat)
mod1 <- svyglm(Yl ~ treatment , design = des1)

summary(mod1)
xtable(blocL[,-6])



for (i in 1:dim(blocL)[1]){
  if (round(bloc.post[i,6],2)==1.00) blocL[i,6] <- "\\checkmark"
}



sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table A6")
print("######")
print("######")
summary(mod1)
blocL
sink()


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Table A8
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

data.st <- read.dta("Data/all data for CH EUvotes _big.dta")
vdata <- data.st
vdata <- cbind(vdata$yes_vote_p,vdata$treated, vdata$kuenstl,vdata$gras,vdata$gebuesch,vdata$baum,vdata$vlose,vdata$wasser,vdata$alti_mun_km,
               vdata$rainfall_surf,vdata$rainfall_1000,vdata$rainfall_1000_2,vdata$Punktflaeche, vdata$SPS_vshare, vdata$CVP_vshare, vdata$GPS_vshare, vdata$FDP_vshare, 
               vdata$SVP_vshare, vdata$vote, vdata$steep_percent, vdata$damageLOW, vdata$damageMIDDLE, vdata$damageHIGH, vdata$elapsedTime, vdata$Gde_typen22, 
               vdata$Kantonsnummer, vdata$Area_50years_km2, vdata$gdenr, vdata$Steuerb_Einkommen, vdata$Steuerpflichtige,vdata$turnouta)
vdata <- data.frame(vdata)
colnames(vdata) <- c("yes_p_per","treated", "kuenstl", "gras", "gebuesch", "baum", "vlose", "wasser", "alti_mun_km", "rainfall_surf", "rainfall_1000", 
                     "rainfall_1000_2", "Punktflaeche", "SPS_vshare", "CVP_vshare", "GPS_vshare", "FDP_vshare", "SVP_vshare", "vote", "steepness",  "damage1", "damage2","damage3","elapsedTime",
                     "Gde_typen22", "cantonnr", "Area_50years_km2", "gdenr", "income", "TaxPopulation", "Turnout")


vdata$kuenst.sh <- 100* vdata$kuenstl/vdata$Punktflaeche
vdata$wasser.sh <- 100* vdata$wasser/vdata$Punktflaeche
vdata$gras.sh <- 100* vdata$gras/vdata$Punktflaeche
vdata$gebuesch.sh <- 100* vdata$gebuesch/vdata$Punktflaeche
vdata$baum.sh <- 100* vdata$baum/vdata$Punktflaeche
vdata$vlose.sh <- 100* vdata$vlose/vdata$Punktflaeche
vdata$alti_mun_m <- 1000* vdata$alti_mun_km
vdata$alti_mun_m2 <- vdata$alti_mun_m^2
vdata$rainfall_surf_1000 <- 1000*vdata$rainfall_surf
vdata$steep_per <- 100 * vdata$steepness

mod1 <- lm(yes_p_per ~  treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare + SVP_vshare +factor(vote) + factor(gdenr),
           data=vdata)
mod2 <- lm(yes_p_per ~ treated +  rainfall_1000  + vlose.sh + wasser.sh + gras.sh + 
             kuenst.sh  +factor(vote) + factor(gdenr), data=vdata)
mod3 <- lm(yes_p_per ~  treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare+ SVP_vshare+  rainfall_1000 +
             vlose.sh + wasser.sh + gras.sh + kuenst.sh  + factor(vote)+ factor(gdenr), data=vdata)

screenreg(list(mod1,mod2,mod3), stars=c(0.01,0.05,0.1), omit="factor|(Intercept)", reorder.coef=c(1:6,7:11),
       custom.coef.names = c("Exposed","Green Party %","Scoial Democrats %" , "Christian Democrats %", 
                             "Liberal Democrats %", "Swiss People's Party %","Rainfall", "No vegetation",
                             "Share of Water", "Share of Gras", "Artificial"), digits=2) 

sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table A8")
print("######")
print("######")
screenreg(list(mod1,mod2,mod3), stars=c(0.01,0.05,0.1), omit="factor|(Intercept)", reorder.coef=c(1:6,7:11),
          custom.coef.names = c("Exposed","Green Party %","Scoial Democrats %" , "Christian Democrats %", 
                                "Liberal Democrats %", "Swiss People's Party %","Rainfall", "No vegetation",
                                "Share of Water", "Share of Gras", "Artificial"), digits=2) 
sink()


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Table A9
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

vdata1 <- na.omit(vdata)
vdata1$vote1 <- vdata1$vote==999001
vdata1$vote2 <- vdata1$vote==999002
vdata1$vote3 <- vdata1$vote==999003
vdata1$vote4 <- vdata1$vote==999004
vdata1$vote5 <- vdata1$vote==999005
vdata1$vote6 <- vdata1$vote==999006
vdata1$vote7 <- vdata1$vote==999007
treatment   <- vdata1$treated 
X           <- cbind(vdata1$vlose.sh, vdata1$wasser.sh, vdata1$gras.sh, vdata1$kuenst.sh, vdata1$alti_mun_m, vdata1$rainfall_surf_1000, 
                     vdata1$steepness, vdata1$SPS_vshare ,vdata1$CVP_vshare,vdata1$GPS_vshare,vdata1$FDP_vshare,vdata1$SVP_vshare,
                     vdata1$vote1, vdata1$vote2, vdata1$vote3, vdata1$vote4, vdata1$vote5, vdata1$vote6,  vdata1$Area_50years_km2)
colnames(X) <- c("no vegetation" ,"wasser", "gras","kuenstl",  "altitude", "rainfall", "steepness", "SPS", "CVP", "GPS", "FDP", "SVP", 
                 "vote1", "vote2", "vote3", "vote4", "vote5", "vote6", "Area_50years_km2")

eb.out <- ebalance(Treatment=treatment, X=X)
weight.L <- rep(1,dim(X)[1])
weight.L[treatment==0] <- eb.out$w
Yl  <- vdata1$yes_p_per
mb.ent.before <- MatchBalance(treatment~ vdata1$Area_50years_km2+vdata1$vlose.sh+ vdata1$wasser.sh+ vdata1$gras.sh  + vdata1$kuenst.sh + 
                                vdata1$alti_mun_m + vdata1$rainfall_surf_1000 + vdata1$steepness + vdata1$SPS_vshare + vdata1$CVP_vshare + 
                                vdata1$GPS_vshare + vdata1$FDP_vshare + vdata1$SVP_vshare, 
                              weights=rep(1,length(weight.L)), 
                              ks = FALSE)
mb.ent.after <- MatchBalance(treatment~ vdata1$Area_50years_km2+vdata1$vlose.sh+ vdata1$wasser.sh+ vdata1$gras.sh  + vdata1$kuenst.sh  + 
                               vdata1$alti_mun_m + vdata1$rainfall_surf_1000 + vdata1$steepness + vdata1$SPS_vshare + vdata1$CVP_vshare + 
                               vdata1$GPS_vshare + vdata1$FDP_vshare + vdata1$SVP_vshare, 
                             weights=weight.L, 
                             ks = FALSE)
bloc.pre <- baltest.collect(matchbal.out=mb.ent.before, var.names=c( "Flooding Risk","Surface: % No vegetation", "Surface: % Water",
                                                                     "Surface: % Gras","Surface: % Artificial","Altitude (in m)", 
                                                                     "Rainfall (per sqkm)", "Steepness in %", "Social Democrats (%)", 
                                                                     "Christian Democrats (%)", "Greens (%)", "Liberals (%)", 
                                                                     "Swiss People's Party (%)"), 
                            after=FALSE)
bloc.post <- baltest.collect(matchbal.out=mb.ent.after, var.names=c("Flooding Risk","Surface: % No vegetation", "Surface: % Water", 
                                                                    "Surface: % Gras","Surface: % Artificial","Altitude (in m)", 
                                                                    "Rainfall (per sqkm)", "Steepness in %", "Social Democrats (%)", 
                                                                    "Christian Democrats (%)", "Greens (%)", "Liberals (%)", 
                                                                    "Swiss People's Party (%)"), 
                             after=FALSE)
rowN <- rownames(bloc.post)
blocL <- matrix(NA,length(rowN),6)
rownames(blocL) <- rowN
blocL[,c(1:3)] <- round(bloc.pre[,c(1,2,6)],2)
blocL[,c(4:5)] <- round(bloc.post[,c(1,2)],2)
dat <- data.frame(Yl,vdata1$cantonnr, weight.L,treatment,vdata1$SPS_vshare,vdata1$GPS_vshare,vdata1$CVP_vshare,
                  vdata1$FDP_vshare,vdata1$SVP_vshare,vdata1$elapsedTime,vdata1$gdenr, vdata1$vote)
des1 <- svydesign(id=~1,weights=~weight.L, data= dat)
mod1 <- svyglm(Yl ~ treatment , design = des1)




for (i in 1:dim(blocL)[1]){
  if (round(bloc.post[i,6],2)==1.00) blocL[i,6] <- "\\checkmark"
}



sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table A9")
print("######")
print("######")
summary(mod1)
blocL
sink()

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Table A10
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

mun.distances <- readRDS("Data/mun_distances.RData")
data.st <- read.dta("Data/all data for CH.dta")
data.st <- data.st[data.st$vote!=226 & data.st$vote!=227 & data.st$vote!=251,]
vdata <- data.st
vdata <- cbind(vdata$yes_vote_meanD,vdata$yes_p_per,vdata$treated, vdata$kuenstl,vdata$gras,vdata$gebuesch,vdata$baum,vdata$vlose,vdata$wasser,vdata$alti_mun_km,
               vdata$rainfall_surf,vdata$rainfall_1000,vdata$rainfall_1000_2,vdata$Punktflaeche, vdata$SPS_vshare, vdata$CVP_vshare, vdata$GPS_vshare, vdata$FDP_vshare, 
               vdata$SVP_vshare, vdata$vote, vdata$steep_percent, vdata$voteL1, vdata$voteL2, vdata$voteL3, vdata$voteL4, vdata$voteL7, vdata$voteL8, vdata$voteL9,
               vdata$voteL10,vdata$voteL11, vdata$damageLOW, vdata$damageMIDDLE, vdata$damageHIGH, vdata$elapsedTime, vdata$small_left_vshare, vdata$Gde_typen22, 
               vdata$Kantonsnummer, vdata$Area_50years_km2, vdata$gdenr, vdata$Tert_share, vdata$Tert_sekun_share, vdata$Steuerb_Einkommen, vdata$Steuerpflichtige, vdata$x_coordinateofcenter,
               vdata$y_coordinateofcenter)
vdata <- data.frame(vdata)
colnames(vdata) <- c("yes_vote_meanD","yes_p_per","treated", "kuenstl", "gras", "gebuesch", "baum", "vlose", "wasser", "alti_mun_km", "rainfall_surf", "rainfall_1000", 
                     "rainfall_1000_2", "Punktflaeche", "SPS_vshare", "CVP_vshare", "GPS_vshare", "FDP_vshare", "SVP_vshare", "vote", "steepness", "voteL1", "voteL2", 
                     "voteL3", "voteL4", "voteL7", "voteL8", "voteL9", "voteL10", "voteL11", "damage1", "damage2","damage3","elapsedTime", "small_left_vshare",
                     "Gde_typen22", "cantonnr", "Area_50years_km2", "gdenr", "Tertiar", "TertiarSekund", "income", "TaxPopulation", "x_coord", "y_coord")#, "Turnout")

vdata$kuenst.sh <- 100* vdata$kuenstl/vdata$Punktflaeche
vdata$wasser.sh <- 100* vdata$wasser/vdata$Punktflaeche
vdata$gras.sh <- 100* vdata$gras/vdata$Punktflaeche
vdata$gebuesch.sh <- 100* vdata$gebuesch/vdata$Punktflaeche
vdata$baum.sh <- 100* vdata$baum/vdata$Punktflaeche
vdata$vlose.sh <- 100* vdata$vlose/vdata$Punktflaeche
vdata$alti_mun_m <- 1000* vdata$alti_mun_km
vdata$alti_mun_m2 <- vdata$alti_mun_m^2
vdata$rainfall_surf_1000 <- 1000*vdata$rainfall_surf
vdata$steep_per <- 100 * vdata$steepness
GDENR <- unique(vdata$gdenr)
'%!in%' <- function(x,y)!('%in%'(x,y))
vdata$treated2 <- NA
vdata$treated4 <- NA
vdata$treated8 <- NA
for (i in unique(vdata$vote)){
  mun.list1 <- list()
  mun.list2 <- list()
  mun.list3 <- list()
  AA <- vdata[vdata$vote==i,]
  treated.mun <- AA$gdenr[AA$treated==1]
  kicker <- which(colnames(mun.distances) %!in%  AA$gdenr )
  mun.distances1 <- mun.distances[-kicker,-kicker]
  #mun.distances1 <- mun.distances
  for (j in treated.mun){
    cL <- which(colnames(mun.distances1)==paste(j))
    treat.distance.j <- mun.distances1[,cL]
    treat.distance.j.bin1 <- rep(0,length(treat.distance.j))
    treat.distance.j.bin2 <- rep(0,length(treat.distance.j))
    treat.distance.j.bin3 <- rep(0,length(treat.distance.j))
    treat.distance.j.bin1[which(treat.distance.j < 2000)] <- 1
    treat.distance.j.bin2[which(treat.distance.j < 4000)] <- 1
    treat.distance.j.bin3[which(treat.distance.j < 8000)] <- 1
    mun.list1[[which(j==treated.mun)]] <- treat.distance.j.bin1
    mun.list2[[which(j==treated.mun)]] <- treat.distance.j.bin2
    mun.list3[[which(j==treated.mun)]] <- treat.distance.j.bin3
  }
  TT1 <- rowSums(matrix(unlist(mun.list1),nrow=nrow(mun.distances1), ncol=length(treated.mun)))
  TT2 <- rowSums(matrix(unlist(mun.list2),nrow=nrow(mun.distances1), ncol=length(treated.mun)))
  TT3 <- rowSums(matrix(unlist(mun.list3),nrow=nrow(mun.distances1), ncol=length(treated.mun)))
  vdata$treated2[vdata$vote==i] <- as.numeric(TT1>0)
  vdata$treated4[vdata$vote==i] <- as.numeric(TT2>0)
  vdata$treated8[vdata$vote==i] <- as.numeric(TT3>0)
  print(i)
}

mod31 <- lm(yes_p_per ~treated2 +  treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare+ SVP_vshare+  rainfall_1000 +
              vlose.sh + wasser.sh + gras.sh + kuenst.sh  + factor(vote) + factor(gdenr), data=vdata)
mod32 <- lm(yes_p_per ~treated4 +  treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare+ SVP_vshare+  rainfall_1000 +
              vlose.sh + wasser.sh + gras.sh + kuenst.sh  + factor(vote) + factor(gdenr), data=vdata)
mod33 <- lm(yes_p_per ~treated8 +  treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare+ SVP_vshare+  rainfall_1000 +
              vlose.sh + wasser.sh + gras.sh + kuenst.sh  + factor(vote) + factor(gdenr), data=vdata)
mod34 <- lm(yes_p_per ~treated2 +treated4 +treated8 +  treated + GPS_vshare + SPS_vshare  + CVP_vshare + FDP_vshare+ SVP_vshare+  rainfall_1000 +
              vlose.sh + wasser.sh + gras.sh + kuenst.sh  + factor(vote) + factor(gdenr), data=vdata)

screenreg(list(mod31,mod32,mod33,mod34), stars=c(0.01,0.05,0.1), omit="factor|(Intercept)", 
       reorder.coef=c(2,1,13:14,3:7,8:12), 
       custom.coef.names = c("Close to Exposure (2 km)","Directly Exposed",
                             "Green Party %","Social Democrats %" , "Christian Democrats %", 
                             "Liberal Democrats %", "Swiss People's Party %","Rainfall", 
                             "No vegetation","Share of Water", "Share of Gras", 
                             "Artificial","Close to Exposure (4 km)","Close to Exposure (8 km)"), digits=2) 

sink("Replication of Tables A.txt", append = TRUE)
print("######")
print("######")
print("Table A10")
print("######")
print("######")
screenreg(list(mod31,mod32,mod33,mod34), stars=c(0.01,0.05,0.1), omit="factor|(Intercept)", 
          reorder.coef=c(2,1,13:14,3:7,8:12), 
          custom.coef.names = c("Close to Exposure (2 km)","Directly Exposed",
                                "Green Party %","Social Democrats %" , "Christian Democrats %", 
                                "Liberal Democrats %", "Swiss People's Party %","Rainfall", 
                                "No vegetation","Share of Water", "Share of Gras", 
                                "Artificial","Close to Exposure (4 km)","Close to Exposure (8 km)"), digits=2) 
sink()





# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Figure 1
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #


gde_15 <- readOGR("geodata/", layer = "gde-1-1-15")
crs(gde_15) <- "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"
gde_15l <- readOGR(dsn="geodata", layer="VEC200_HOHEITSGEBIET_LV95")
gde_15l <- spTransform(gde_15l, CRS("+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"))
ch <- gde_15l[gde_15l@data$ICC=="CH", ]
ch_geom <- fortify(ch, region="ICC")
seen_geom <- readRDS("geodata/seen_geom.Rds")
ch_geom <- readRDS("geodata/ch_geom.Rds")
# read in event data, coded as "quakes"
quakes <- foreign::read.dta("geodata/Disaster_events_dates_keeperL.dta")%>%
  dplyr::select(J, K, keeperL)%>%
  rename(long=J, lat=K)
quakes0 <- quakes[quakes$keeperL==0,1:2]
quakes1 <- quakes[quakes$keeperL==1,1:2]
map_data <- readRDS("geodata/map_data.Rds")
#read in relief ----
relief <- raster("geodata/02-relief-georef-clipped-resampled.tif")
relief_spdf <- as(relief, "SpatialPixelsDataFrame")
relief <- as.data.frame(relief_spdf) %>% 
  rename(value = `X02.relief.georef.clipped.resampled`)
# actual map
colls <- c("violetred4", "firebrick", "firebrick4", "violetred", "darkred", "darkmagenta", "darkorchid4", "blueviolet", "deeppink1")
col0 <- "white"
col1 <- "firebrick"
i <- 1
ggplot()+
  geom_raster(data = relief, aes(x = x, y = y, alpha=value/100))+
  geom_polygon(data = map_data, aes(x=long, y=lat, group=group), fill="white")+ 
  geom_path(data = ch_geom, aes(x=long, y=lat, group=group),color="grey80", size=.3)+
  geom_polygon(data=seen_geom, aes(x=long, y=lat, group=group), fill="lightblue")+  
  geom_path(data=seen_geom, aes(x=long, y=lat, group=group), color="grey80", size=.3)+
  geom_point(data=quakes0, aes(x=long, y=lat), color=col0, alpha=.5, size=0.4, pch=4)+
  geom_point(data=quakes1, aes(x=long, y=lat), color=col1, alpha=.8, size=0.8)+
  coord_equal()+
  guides(alpha=FALSE)+
  theme_bw()+
  labs(y="", x="")+ 
  theme(axis.text=element_blank(),
        axis.ticks = element_blank(),
        panel.grid = element_line(size=.2))+
  ggsave(filename="Relief_map.pdf", height=6, width=9.708)








